#MT Analysis of Survey Disposition
#PORES AAPOR
#December 2020

rm(list=ls())
#Set to the folder you've downloaded the material to
setwd()
load(file="DispositionDataForAnalysis.Rdata")
library(lfe)
library(stargazer)
library(survey)


###########
#Appendix A : Summary statistics table
##########


keep <- c("coop", "contact", "dem","rep", "ind",
          "female", "AgeUnder30", "Age3039", "Age4049", "Age5059",
          "Age4049", "Age5059", "Age6074", "Age75p",
          "white", "black", "hispanic", "other.race", "turnout20",
          "dem.perc16","urban", "rural","suburban","cases" )

dat$dem[dat$likely.party=="D"] <- 1
dat$dem[dat$likely.party %in% c("R","I")] <- 0
dat$rep[dat$likely.party=="R"] <- 1
dat$rep[dat$likely.party %in% c("D","I")] <- 0
dat$ind[dat$likely.party=="I"] <- 1
dat$ind[dat$likely.party %in% c("D","R")] <- 0

dat$rural[dat$density=="R"] <- 1
dat$rural[dat$density %in% c("U","S")] <- 0
dat$urban[dat$density=="U"] <- 1
dat$urban[dat$density %in% c("R","S")] <- 0
dat$suburban[dat$density=="S"] <- 1
dat$suburban[dat$density %in% c("U","R")] <- 0



keep %in% names(dat)
sum <- dat[keep]

stargazer(sum, summary.stat =c("n", "mean", "sd", "min", "max"))

#Partisan raw numbers in the disposition file

state <- unique(dat$state)
partisan.distribution <- matrix(NA, nrow=length(state), ncol=6)

partisan.distribution[,1] <- state
partisan.distribution

for(i in 1:length(state)){
  partisan.distribution[i,2] <- nrow(dat[dat$likely.party=="D" & dat$state==state[i] & !is.na(dat$likely.party),])
  partisan.distribution[i,3] <- nrow(dat[dat$likely.party=="I" & dat$state==state[i],])
  partisan.distribution[i,4] <- nrow(dat[dat$likely.party=="R" & dat$state==state[i],])
}


partisan.distribution[,5] <- c("Imputed", "Registration", "Registration", "Registration", "Registration", "Imputed", "Imputed", "Registration", "Registration", "Registration", "Imputed", "Imputed")
partisan.distribution[,6] <- c("No", "Yes","Yes", "No", "Yes", "Yes", "Yes", "No", "No", "Yes", "No", "Yes")

library(xtable)
print(xtable(partisan.distribution), include.rownames=F)



#Partisan raw numbers in the exit polls file

state <- unique(ep$state)

partisan.distribution <- matrix(NA, nrow=length(state), ncol=4)

partisan.distribution[,1] <- state
partisan.distribution

for(i in 1:length(state)){
  partisan.distribution[i,2] <- nrow(ep[ep$likely.party=="D" & ep$state==state[i],])
  partisan.distribution[i,3] <- nrow(ep[ep$likely.party=="I" & ep$state==state[i],])
  partisan.distribution[i,4] <- nrow(ep[ep$likely.party=="R" & ep$state==state[i],])
}



library(xtable)
print(xtable(partisan.distribution), include.rownames=F)


#Figure of distribution of partisanship in exit poll
# raw, weighted with original and weighted with re-weight
state <- unique(ep$state)
state <- rev(state)

raw <- matrix(NA, nrow=length(state), ncol=3)
w1  <- matrix(NA, nrow=length(state), ncol=3)
w2  <- matrix(NA, nrow=length(state), ncol=3)

for(j in 1:length(state)){
  raw[j,] <- prop.table(table(ep$likely.party[ep$state==state[j]]))*100

  #Generate state disposition data
  d <- dat[dat$state==state[j],]
  #Generate state exit poll data
  e <- ep[ep$state==state[j],]
  #Generate additional id, to play well with survey package
  e$id2 <- seq(1, nrow(e))

  #Use logit to get cooperation rates by party
  #(This returns the exact same result as calculating the means,
  #but using logit allows this to be scaled to more variables.)
  m <- glm(coop ~  likely.party, data=d, family = "binomial")

  #Make reduced dataset for prediction
  keep <- c("id2","likely.party")
  e.p<- e[keep]
  e.p <-e.p[complete.cases(e.p),]
  #Predict cooperation rate in the exit poll using logit result
  e.p$predicted.coop <- predict(m, type="response", newdata=e.p)

  #Invert predicted cooperation to get a cooperation weight
  e.p$coop.weight <- 1/(e.p$predicted.coop)

  #Merge this back into exitpoll data
  keep <- c("id2","coop.weight")
  e.p <- e.p[keep]
  e <- merge(e,e.p, by="id2")

  #Generate weighting objects
  original.weights <- svydesign(id=~id2, weights=~weight, data=e)
  e$reweight.coop <- e$weight*e$coop.weight
  coop.weights <- svydesign(id=~id2, weights=~reweight.coop, data=e)

  w1[j,] <- prop.table(svytable(~likely.party, design=original.weights))*100
  w2[j,] <- prop.table(svytable(~likely.party, design=coop.weights))*100
}

ypos <- seq(1,length(state))
ypos2 <- ypos-.25
ypos3 <- ypos - .5

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureA1.pdf", height=6, width=6)
plot(raw[,1], ypos, col="dodgerblue", pch=16, type="n", xlim=c(0,100), ylim=c(0,7),
     axes=F, ylab="", xlab="Percent")
abline(v=50, lty=2)

points(raw[,1], ypos, col="dodgerblue", pch=16)
points(raw[,2], ypos, col="purple", pch=16)
points(raw[,3], ypos, col="firebrick", pch=16)


points(w1[,1], ypos2, col="dodgerblue", pch=17)
points(w1[,2], ypos2, col="purple", pch=17)
points(w1[,3], ypos2, col="firebrick", pch=17)

points(w2[,1], ypos3, col="dodgerblue", pch=18)
points(w2[,2], ypos3, col="purple", pch=18)
points(w2[,3], ypos3, col="firebrick", pch=18)
abline(h=ypos-.75, lty=2)
axis(side=1)
axis(side=2, at=seq(1,7,1), labels=state, las=2)
legend("topright", c("Raw in Exit Poll", "Original Weights", "Coop Re-Weight"), pch=c(16,17,18), bg="white")
dev.off()




###########
#Appendix B : Full Results and alternative specifications
##########


#Cooperation

#Main specification
m1 <- felm(coop ~  likely.party + log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic   | state | 0 | state, data=dat)


#Simple specification
m2 <- felm(coop ~  likely.party   | state | 0 | state, data=dat)

#Individual Demos Only
m3 <- felm(coop ~  likely.party +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic      | state | 0 | state, data=dat)


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")


#Full Results for Appendix
stargazer(m2, m3,m1,  style="apsr",
          dep.var.labels = "P(Cooperate)", type="text")
stargazer(m2, m3,m1, style="apsr",
          dep.var.labels = "P(Cooperate)",
          covariate.labels = labels,
          column.labels = c("Party Only", "Individual Controls", "Full Specification"),
          add.lines= list(c("State FE", "Yes","Yes","Yes")),
          digits.extra = 1
          )



#Coontact

#Main specification
m1 <- felm(contact ~  likely.party + log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic   | state | 0 | state, data=dat)


#Simple specification
m2 <- felm(contact ~  likely.party   | state | 0 | state, data=dat)

#Individual Demos Only
m3 <- felm(contact ~  likely.party +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic      | state | 0 | state, data=dat)


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")


#Full Results for Appendix
stargazer(m2, m3,m1,  style="apsr",
          dep.var.labels = "P(Contact)", type="text")
stargazer(m2, m3,m1, style="apsr",
          dep.var.labels = "P(Contact)",
          covariate.labels = labels,
          column.labels = c("Party Only", "Individual Controls", "Full Specification"),
          add.lines= list(c("State FE", "Yes","Yes","Yes")),
          digits.extra = 1
)


#Overall Response

dat$response <- NA
dat$response[dat$dispo.cleaned==1] <- 1
dat$response[dat$dispo.cleaned %in% c(2,3,4)] <- 0
table(dat$response)

#Main specification
m1 <- felm(response ~  likely.party + log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic   | state | 0 | state, data=dat)


#Simple specification
m2 <- felm(response ~  likely.party   | state | 0 | state, data=dat)

#Individual Demos Only
m3 <- felm(response ~  likely.party +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic      | state | 0 | state, data=dat)


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")


#Full Results for Appendix
stargazer(m2, m3,m1,  style="apsr",
          dep.var.labels = "P(Response)", type="text")
stargazer(m2, m3,m1, style="apsr",
          dep.var.labels = "P(Response)",
          covariate.labels = labels,
          column.labels = c("Party Only", "Individual Controls", "Full Specification"),
          add.lines= list(c("State FE", "Yes","Yes","Yes")),
          digits.extra = 1
)



#Overall Response, excluding DQ

dat$response <- NA
dat$response[dat$dispo.cleaned==1] <- 1
dat$response[dat$dispo.cleaned %in% c(2,3)] <- 0
table(dat$response)

#Main specification
m1 <- felm(response ~  likely.party + log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic   | state | 0 | state, data=dat)


#Simple specification
m2 <- felm(response ~  likely.party   | state | 0 | state, data=dat)

#Individual Demos Only
m3 <- felm(response ~  likely.party +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic      | state | 0 | state, data=dat)


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")


#Full Results for Appendix
stargazer(m2, m3,m1,  style="apsr",
          dep.var.labels = "P(Response)", type="text")
stargazer(m2, m3,m1, style="apsr",
          dep.var.labels = "P(Response|Non-DQ)",
          covariate.labels = labels,
          column.labels = c("Party Only", "Individual Controls", "Full Specification"),
          add.lines= list(c("State FE", "Yes","Yes","Yes")),
          digits.extra = 1
)





###########
#Appendix C : State-by State Results
##########

#Cooperate Table

#Creating table of every state regression output
state <- sort(unique(dat$state[!is.na(dat$likely.party)]))

m1 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[1],])
m2 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[2],])
m3 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[3],])
m4 <- felm(coop ~likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[4],])
m5 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[5],])
m6 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[6],])
m7 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[7],])
m8 <- felm(coop ~likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[8],])
m9 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[9],])
m10 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic , data=dat[dat$state==state[10],])
m11 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic , data=dat[dat$state==state[11],])
m12 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic, data=dat[dat$state==state[12],])


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")

stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,
          covariate.labels = labels,
          column.labels = state,
          dep.var.labels = "P(Cooperate)",
          model.numbers = F)
#Contact Table

state <- sort(unique(dat$state[!is.na(dat$likely.party)]))

m1 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[1],])
m2 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[2],])
m3 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[3],])
m4 <- felm(contact ~likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[4],])
m5 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[5],])
m6 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[6],])
m7 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic, data=dat[dat$state==state[7],])
m8 <- felm(contact ~likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[8],])
m9 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic , data=dat[dat$state==state[9],])
m10 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic , data=dat[dat$state==state[10],])
m11 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic , data=dat[dat$state==state[11],])
m12 <- felm(contact ~ likely.party + log(cases) + dem.perc16 + density +
              female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic, data=dat[dat$state==state[12],])


labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")

stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,
          covariate.labels = labels,
          column.labels = state,
          dep.var.labels = "P(Contact)",
          model.numbers = F)


#Contact Graph
state <- sort(unique(dat$state[!is.na(dat$likely.party)]))
rep.contact <- NA
ind.contact <- NA
rep.contact.se <- NA
ind.contact.se <- NA

for(i in 1:length(state)){
  m1 <- felm(contact ~ likely.party + dem.perc16 + log(cases) +
               female +AgeUnder30  +Age3039 + Age4049 + Age5059 + Age6074 +
               white + black + hispanic +density, data=dat[dat$state==state[i],])

  rep.contact[i] <- coef(m1)["likely.partyR"]
  ind.contact[i]  <- coef(m1)["likely.partyI"]
  rep.contact.se[i]  <- sqrt(vcov(m1)["likely.partyR","likely.partyR"])
  ind.contact.se[i]  <-  sqrt(vcov(m1)["likely.partyI","likely.partyI"])
}


ypos1 <- seq(12,1,-1)
ypos2 <- seq(11.5,0.5,-1)
ypos3 <- seq(11.75,.75,-1)

#Contact
pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureC1.pdf")
plot(rep.contact,ypos1, xlim=c(-0.20,0.20), ylim=c(0,13), axes=F, pch=16, col="firebrick",
     ylab="", xlab="Effect on P(Contact) vs. Democrat")
segments(rep.contact,ypos1, rep.contact +rep.contact.se*1.96,ypos1, col="firebrick")
segments(rep.contact,ypos1, rep.contact -rep.contact.se*1.96,ypos1, col="firebrick")


points(ind.contact, ypos2,pch=16, col="purple")
segments(ind.contact,ypos2, ind.contact +ind.contact.se*1.96,ypos2, col="purple")
segments(ind.contact,ypos2, ind.contact -ind.contact.se*1.96,ypos2, col="purple")


axis(side=1, at=seq(-0.20, 0.20,0.02))
axis(side=2, at=seq(11.75,.75,-1), labels=state, las=2)
abline(h=seq(1.25,11.25,1),lty=2, col="gray80")
abline(v=0, lty=3)
legend("topright", c("Republican","Independent/Other"), pch=c(16,16), col=c("firebrick","purple"))
dev.off()




###########
#Appendix D : Comparison of Registration and Imputation Party States
##########

#Does imputing party affect things?


m1 <- felm(coop ~  likely.party +  log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic + likely.party| state | 0 | state, data=dat[dat$imputed.party==0,])

m2 <- felm(coop ~ likely.party + log(cases) + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic | state | 0 | state, data=dat[dat$imputed.party==1,])



stargazer(m1,m2, type="text")

labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic")

stargazer(m1,m2, style="apsr",
          dep.var.labels = "P(Cooperate)",
          column.labels = c("Party Registration States","Imputed Party States"),
          covariate.labels = labels)


###########
#Appendix E : Comparison of Original and Cooperation-Adjusted Weights
##########

#States
state <- sort(unique(ep$state))
state <- rev(state)

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureE1.pdf", height=11, width=8.5)
par(mfrow=c(4,2))
for(j in 1:length(state)){
  #Generate state disposition data
  d <- dat[dat$state==state[j],]
  #Generate state exit poll data
  e <- ep[ep$state==state[j],]
  #Generate additional id, to play well with survey package
  e$id2 <- seq(1, nrow(e))
  #Use logit to get cooperation rates by party
  #(This returns the exact same result as calculating the means,
  #but using logit allows this to be scaled to more variables.)
  m <- glm(coop ~  likely.party, data=d, family = "binomial")

  #Make reduced dataset for prediction
  keep <- c("id2","likely.party")
  e.p<- e[keep]
  e.p <-e.p[complete.cases(e.p),]
  #Predict cooperation rate in the exit poll using logit result
  e.p$predicted.coop <- predict(m, type="response", newdata=e.p)

  #Invert predicted cooperation to get a cooperation weight
  e.p$coop.weight <- 1/(e.p$predicted.coop)

  #Merge this back into exitpoll data
  keep <- c("id2","coop.weight")
  e.p <- e.p[keep]
  e <- merge(e,e.p, by="id2")

  #Generate weighting objects
  e$reweight.coop <- e$weight*e$coop.weight

  #plot
  plot(e$weight[e$likely.party=="I"], e$reweight.coop[e$likely.party=="I"], xlab="Original Weight",
       ylab="Cooperation Re-Weight", pch=1,
       main=print(state[j]), col="purple")
  points(e$weight[e$likely.party=="D"], e$reweight.coop[e$likely.party=="D"],  pch=1,
         col="dodgerblue")
  points(e$weight[e$likely.party=="R"], e$reweight.coop[e$likely.party=="R"],  pch=1,
         col="firebrick")
}
dev.off()


###########
#Appendix F : Cooperation Adjustment with Full Demographics
##########

#States
state <- sort(unique(ep$state))
state <- rev(state)
margin.weighted1 <- rep(NA, length(state))
margin.weighted2<- rep(NA, length(state))
trump.actual <- rep(NA, length(state))
biden.actual <- rep(NA, length(state))

for(j in 1:length(state)){
  trump.actual[j] <- el$rep_percent[el$stateid==state[j]]
  biden.actual[j] <- el$dem_percent[el$stateid==state[j]]
}

trump.actual <- as.numeric(gsub("[\\%]","",trump.actual))
biden.actual <- as.numeric(gsub("[\\%]","",biden.actual))
margin.real <- biden.actual - trump.actual
survey.n <- rep(NA, length(state))

#Reweighting of polls, by state

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureF2.pdf", height=11, width=8.5)
par(mfrow=c(4,2))
for(j in 1:length(state)){
  #Generate state disposition data
  d <- dat[dat$state==state[j],]
  #Generate state exit poll data
  e <- ep[ep$state==state[j],]
  #Generate additional id, to play well with survey package
  e$id2 <- seq(1, nrow(e))
  #Save n of survey in state
  survey.n[j] <- nrow(e)

  #Use logit to build model for probability of cooperation.
  #This is step that differs from above
  m <- glm(coop ~  likely.party +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             white + black + hispanic , data=d, family = "binomial")

  summary(m)

  #Make reduced dataset for prediction

  keep <- c("id2","likely.party","female","Age3039","Age4049","Age5059","Age6074",
            "Age75p","white","black","hispanic")
  e.p<- e[keep]
  e.p <-e.p[complete.cases(e.p),]

  #Predict cooperation rate in the exit poll using logit result
  e.p$predicted.coop <- predict(m, type="response", newdata=e.p)

  #Invert predicted cooperation to get a cooperation weight
  e.p$coop.weight <- 1/(e.p$predicted.coop)

  #Merge this back into exitpoll data
  keep <- c("id2","coop.weight")
  e.p <- e.p[keep]
  e <- merge(e,e.p, by="id2")

  #Generate weighting objects
  original.weights <- svydesign(id=~id2, weights=~weight, data=e)
  e$reweight.coop <- e$weight*e$coop.weight
  coop.weights <- svydesign(id=~id2, weights=~reweight.coop, data=e)

  #Original Estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[1]
  margin.weighted1[j] <- biden - trump

  #Cooperation adjusted estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=coop.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=coop.weights)))[1]
  margin.weighted2[j] <- biden - trump
  #plot
  plot(e$weight[e$likely.party=="I"], e$reweight.coop[e$likely.party=="I"], xlab="Original Weight",
       ylab="Cooperation Re-Weight", pch=1,
       main=print(state[j]), col="purple")
  points(e$weight[e$likely.party=="D"], e$reweight.coop[e$likely.party=="D"],  pch=1,
         col="dodgerblue")
  points(e$weight[e$likely.party=="R"], e$reweight.coop[e$likely.party=="R"],  pch=1,
         col="firebrick")
}
dev.off()

margin.weighted1 <- round(margin.weighted1*100,2)
margin.weighted2 <- round(margin.weighted2*100,2)
margin.real <- round(margin.real,2)

ypos <- seq(1,length(state))

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureF1.pdf", height=6, width=6)
plot(margin.real,ypos, xlim=c(-30,30), pch=16, xlab="Estimate of Biden - Trump Margin", ylab="",
     axes=F, col="black", type="p")
points(margin.weighted1, ypos, pch=17, col="darkorange")
points(margin.weighted2, ypos, pch=18, col="forestgreen")
segments(margin.weighted1, ypos,margin.weighted2, ypos, col="forestgreen", lwd=2)

abline(v=0, lty=2)
axis(side=1, at=seq(-30,30,10))
axis(side=2, at=ypos, labels=state, las=2)
legend("topright", c("Election Result","Original Survey Estimate", "Coop Reweighted Estimate"),
       pch=c(16,17,18), col=c("black","darkorange","forestgreen"), cex=.75)
dev.off()


improve <- (margin.weighted1 - margin.real) - (margin.weighted2 - margin.real)
weighted.mean(improve, w=1/survey.n)



###########
#Appendix G : Relationship between Likely Party and Partisan Self Identification
##########




#Full JPF
library(kableExtra)

#AZ,CO,IA,MI,MN,PA,WI
table(ep$vote.choice)
ep$biden <- NA
ep$biden[ep$vote.choice=="Biden"]<- 1
ep$biden[ep$vote.choice=="Trump"]<- 0

likely.party <- c("D","I","R")
pid.lab <- c("Dem","Ind","Rep")

#Overall

round(prop.table(table(ep$pid.lab, ep$likely.party),2)*100,1)


vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i]],na.rm=T)
  }
}
round(vote.matrix*100,2)

mean(ep$biden[ep$likely.party=="D"],na.rm=T)
mean(ep$biden[ep$likely.party=="I"],na.rm=T)
mean(ep$biden[ep$likely.party=="R"],na.rm=T)

#AZ
round(prop.table(table(ep$pid.lab[ep$state=="AZ"], ep$likely.party[ep$state=="AZ"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="AZ"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)


#CO
round(prop.table(table(ep$pid.lab[ep$state=="CO"], ep$likely.party[ep$state=="CO"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="CO"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)

#IA
round(prop.table(table(ep$pid.lab[ep$state=="IA"], ep$likely.party[ep$state=="IA"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="IA"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)

#MI
round(prop.table(table(ep$pid.lab[ep$state=="MI"], ep$likely.party[ep$state=="MI"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="MI"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)

#MN
round(prop.table(table(ep$pid.lab[ep$state=="MN"], ep$likely.party[ep$state=="MN"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="MN"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)


#PA
round(prop.table(table(ep$pid.lab[ep$state=="PA"], ep$likely.party[ep$state=="PA"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="PA"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)

#WI
round(prop.table(table(ep$pid.lab[ep$state=="WI"], ep$likely.party[ep$state=="WI"]),2)*100,2)

vote.matrix <- matrix(NA, nrow=3, ncol=3)
for(i in 1:3){
  for(j in 1:3){
    vote.matrix[i,j] <- mean(ep$biden[ep$likely.party==likely.party[j] & ep$pid.lab==pid.lab[i] & ep$state=="WI"] ,na.rm=T)
  }
}
round(vote.matrix*100,2)




###########
#Appendix H : Weighting to Distribution of Partisanship in the Voter File
##########

state <- sort(unique(ep$state))
state <- rev(state)

margin.weighted1 <- rep(NA, length(state))
margin.weighted2<- rep(NA, length(state))
trump.actual <- rep(NA, length(state))
biden.actual <- rep(NA, length(state))

for(j in 1:length(state)){
  trump.actual[j] <- el$rep_percent[el$stateid==state[j]]
  biden.actual[j] <- el$dem_percent[el$stateid==state[j]]
}

trump.actual <- as.numeric(gsub("[\\%]","",trump.actual))
biden.actual <- as.numeric(gsub("[\\%]","",biden.actual))
margin.real <- biden.actual - trump.actual
survey.n <- rep(NA, length(state))

for(j in 1:length(state)){
  e <- ep[ep$state==state[j],]
  e <- e[!is.na(e$weight),]
  e <- e[!is.na(e$pid),]
  survey.n[j] <- nrow(e)
  e$id2 <- seq(1,nrow(e))
  #Generate new weight
  original.weights <- svydesign(id=~id2, weights=~weight, data=e)
  v <- vf.id[vf.id$state==state[j],]
  v$state <- NULL


  party.weights <- postStratify(original.weights, strata = ~likely.party, population =v)

  prop.table(svytable(~likely.party, design=original.weights))
  prop.table(svytable(~likely.party, design=party.weights))


  #Original Estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[1]
  margin.weighted1[j] <- biden - trump

  #Party adjusted estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=party.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=party.weights)))[1]
  margin.weighted2[j] <- biden - trump
}


margin.weighted1 <- round(margin.weighted1*100,2)
margin.weighted2 <- round(margin.weighted2*100,2)
margin.real <- round(margin.real,2)

ypos <- seq(1,length(state))

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureH1.pdf", height=6, width=6)
plot(margin.real,ypos, xlim=c(-30,30), pch=16, xlab="Estimate of Biden - Trump Margin", ylab="",
     axes=F, col="black", type="p")
points(margin.weighted1, ypos, pch=16, col="darkorange")
points(margin.weighted2, ypos, pch=16, col="forestgreen")
segments(margin.weighted1, ypos,margin.weighted2, ypos, col="forestgreen", lwd=2)

abline(v=0, lty=2)
axis(side=1, at=seq(-30,30,10))
axis(side=2, at=ypos, labels=state, las=2)
legend("topright", c("Election Result","Original Survey Estimate", "Party Reweighted Estimate"),
       pch=c(16,16,16), col=c("black","darkorange","forestgreen"), cex=.75)
dev.off()


improve <- (margin.weighted1 - margin.real) - (margin.weighted2 - margin.real)
weighted.mean(improve, w=1/survey.n)
cbind(state,survey.n)
cbind(state,improve)


#Comparison of partisanship across sources
#Create table of disposition likely partisanship versus voter file


state <- c("AZ","CO","IA","MI","MN","PA","WI")
state <- rev(state)

dem.vf <- rep(NA, length(state))
ind.vf <- rep(NA, length(state))
rep.vf <- rep(NA, length(state))
dem.ep <- rep(NA, length(state))
ind.ep <- rep(NA, length(state))
rep.ep <- rep(NA, length(state))
dem.disp <- rep(NA, length(state))
ind.disp <- rep(NA, length(state))
rep.disp <- rep(NA, length(state))

for(i in 1:length(state)){
  dem.vf[i] <- vf.id$perc[vf.id$state==state[i] & vf.id$likely.party=="D"]
  ind.vf[i] <- vf.id$perc[vf.id$state==state[i] & vf.id$likely.party=="I"]
  rep.vf[i] <- vf.id$perc[vf.id$state==state[i] & vf.id$likely.party=="R"]
  dem.ep[i] <- ep.id$D[ep.id$state==state[i]]
  ind.ep[i] <- ep.id$I[ep.id$state==state[i]]
  rep.ep[i] <- ep.id$R[ep.id$state==state[i]]
  dem.disp[i]<- prop.table(table(dat$likely.party[dat$state==state[i]]))[1]*100
  ind.disp[i]<- prop.table(table(dat$likely.party[dat$state==state[i]]))[2]*100
  rep.disp[i]<- prop.table(table(dat$likely.party[dat$state==state[i]]))[3]*100
}


ypos <- seq(1,length(state))
ypos2 <- ypos-.25
ypos3 <- ypos - .5

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureH2.pdf", height=6, width=6)
plot(dem.vf, ypos, col="dodgerblue", pch=16,
     xlim=c(0,100), ylim=c(0,7),xlab="Percent",
     ylab="", axes=F)
points(rep.vf, ypos, col="firebrick", pch=16)
points(ind.vf, ypos, col="purple", pch=16)
points(dem.ep, ypos2, col="dodgerblue", pch=17)
points(rep.ep, ypos2, col="firebrick", pch=17)
points(ind.ep, ypos2, col="purple", pch=17)
points(dem.disp, ypos3, col="dodgerblue", pch=18)
points(rep.disp, ypos3, col="firebrick", pch=18)
points(ind.disp, ypos3, col="purple", pch=18)

abline(h=ypos-.8, lty=2)
axis(side=2, at=ypos-.125, labels=state, las=2)
axis(side=1, at=seq(0,100,20))
legend("topright", c("Voter File","Weight to Outcome", "Disposition File"), pch=c(16,17,18), bg = "white")
dev.off()






###########
#Appendix I : Investigating impact of COVID-19
##########



m1 <- felm(contact ~  log(cases)*likely.party + dem.perc16 + density +
             female +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic + likely.party  | state | 0 | state, data=dat)

m2 <- felm(coop ~  log(cases)*likely.party + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic  | state | 0 | state, data=dat)

stargazer(m1,m2, type="text")
getfe(m1)

labels <- c("ln(COVID Cases)","Independent/Other","Republican","Dem \% Pres. Election 2020",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic","ln(COVID Cases)*Independent/Other","ln(COVID Cases)*Republican")

#Full Results for Appendix
stargazer(m1,m2,  style="apsr",
          dep.var.labels = c("P(Contact)","P(Cooperate)"),
          covariate.labels = labels)



summary(dat$cases)

#Custom margins function
Rmargins <- function(model,xname,xzname,vcov.matrix,levels, ci=.95, latex=F){
  # Generate Marginal Effects
  m_effect <- coef(model)[xname]+coef(model)[xzname]*levels
  #Variances and Covariances
  v_x <- vcov.matrix[xname,xname]
  v_xz <- vcov.matrix[xzname,xzname]
  cov <- vcov.matrix[xname,xzname]
  #Standard Errors
  se <- sqrt(v_x + (levels^2)*v_xz+2*levels*cov)
  #T value for 95%CI
  if (ci==.95){
    t <- qt(0.025,model$df)
  }
  #T value for 90%CI
  if (ci==.9){
    t <- qt(0.05,model$df)
  }
  #Confidence Bounds
  m_upper <- m_effect + se*t
  m_lower <- m_effect - se*t
  # Remove Flotsom and Jetson
  #Printing Table
  table <- cbind(levels,m_effect,se,m_upper,m_lower)
  if (latex==T){
    library(xtable)
    print(xtable(table),include.rownames=FALSE)
  }
  if (latex==F){
    return(table)
  }
}


cases <- seq(1000,50000, 10)

m.i <- Rmargins(m1, "likely.partyI", "log(cases):likely.partyI", vcov(m1), levels=log(cases) )
m.r <- Rmargins(m1, "likely.partyR", "log(cases):likely.partyR", vcov(m1), levels=log(cases) )
cases.sample <- sample(dat$cases, 1000)


pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureI1.pdf", height=6, width=6)
plot(cases, m.i[,2], col="purple", lwd=4, type="n", ylim=c(-.01, .01), ylab="Effect of Likely Partisanship, Relative to Democrat", xlab="County COVID Cases",
     main="P(Contact)")
abline(h=0, lty=2, lwd=2)
points(cases, m.i[,2], type="l", col="purple", lwd=4)
points(cases, m.i[,4], col="purple", lty=2, type="l")
points(cases, m.i[,5], col="purple", lty=2, type="l")

points(cases, m.r[,2], col="firebrick", type="l", lwd=4)
points(cases, m.r[,4], col="firebrick", lty=2, type="l")
points(cases, m.r[,5], col="firebrick", lty=2, type="l")
rug(cases.sample)
legend("topleft", c("Republican", "Independent/Other"), lty=c(1,1), col=c("firebrick","purple"))
dev.off()

m.i <- Rmargins(m2, "likely.partyI", "log(cases):likely.partyI", vcov(m2), levels=log(cases) )
m.r <- Rmargins(m2, "likely.partyR", "log(cases):likely.partyR", vcov(m2), levels=log(cases) )

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER FigureI2.pdf", height=6, width=6)
plot(cases, m.i[,2], col="purple", lwd=4, ylim=c(-.1,0),type="n", ylab="Effect of Likely Partisanship, Relative to Democrat", xlab="County COVID Cases",
     main="P(Cooperate)")
abline(h=0, lty=2, lwd=2)
points(cases, m.i[,2], type="l", col="purple", lwd=4)
points(cases, m.i[,4], col="purple", lty=2, type="l")
points(cases, m.i[,5], col="purple", lty=2, type="l")

points(cases, m.r[,2], col="firebrick", type="l", lwd=4)
points(cases, m.r[,4], col="firebrick", lty=2, type="l")
points(cases, m.r[,5], col="firebrick", lty=2, type="l")
legend("bottomright", c("Republican", "Independent/Other"), lty=c(1,1), col=c("firebrick","purple"))
rug(cases.sample)

dev.off()




###########
#Appendix K : Investigating confounding effects of turnout
##########

#Cooperation

#2016 Turnout
m16 <- felm(coop ~  likely.party + log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic + turnout16   | state | 0 | state, data=dat)

#2020 Turnout
m20 <- felm(coop ~  likely.party + log(cases) + dem.perc16 + density +
              female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic + turnout20   | state | 0 | state, data=dat)


#Summed turnout
dat$turnout.sum <- dat$turnout16 + dat$turnout20

mb <- felm(coop ~  likely.party + log(cases) + dem.perc16 + density +
              female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
              other.race + black + hispanic + turnout.sum   | state | 0 | state, data=dat)



labels <- c("Independent/Other","Republican","ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic", "Voted in 2016 Election", "Voted in 2020 Election", "Summed Turnout")


#Full Results for Appendix
stargazer(m16,m20, mb, style="apsr",
          dep.var.labels = "P(Cooperate)", type="text")
stargazer(m16,m20, mb, style="apsr",
          dep.var.labels = "P(Cooperate)",
          covariate.labels = labels,
          add.lines= list(c("State FE", "Yes","Yes")),
          digits.extra = 1
)





